{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiline TRL"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiline TRL is a two-port VNA calibration utilizing at least two transmission lines with different physical lengths and at least one reflective standard that is identical on both ports. The electrical parameters of the lines don't need to be known, but the transmission lines should have identical construction (same propagation constant and characteristic impedance). The reflect standard reflection coefficient doesn't need to be known exactly, phase needs to be known with 90 degree accuracy.\n",
    "\n",
    "If the measured phase differences of the lines is a multiple of 180 degrees the calibration is singular. The calibration accuracy is worse the closer the line measurement phases are to the singularities, the best accuracy is obtained in the two lines case when the phase difference is 90 degrees. Multiple lines can be used to extend the frequency range where the calibration is accurate."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example demonstrates  how to use `skrf`'s NIST-style Multiline calibration (`NISTMultilineTRL`). First a [simple application](#Simple-Multiline-TRL) is presented, followed by a [full simulation](#Compare-calibrations-with-different-combinations-of-lines) to demonstrate the improvements in calibration accuracy vs the number of lines.  All  data is used in the demonstration is generated by skrf, and the code for this is given [at the end of the example](#Simulation-to-generate-the-input-data)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition, this example showcases the TUG multiline calibration (TUGMultilineTRL), which is an alternative to the commonly used NIST-style multiline method. The TUG method offers better statistical properties. An example also provided comparing both calibration methods through [Monte Carlo analysis](#monte-carlo-analysis---nist-vs-tug)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Simple Multiline TRL"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import skrf\n",
    "from skrf.calibration import NISTMultilineTRL, TUGMultilineTRL, terminate\n",
    "from skrf.media import CPW, Coaxial\n",
    "from skrf.network import two_port_reflect\n",
    "\n",
    "skrf.stylely()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Load data into skrf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load all measurement data into a dictionary\n",
    "data = skrf.io.general.read_all_networks('multiline_trl_data/')\n",
    "\n",
    "# Pull out measurements by name into an ordered list\n",
    "measured_names = ['thru','reflect','linep3mm','line2p3mm']\n",
    "measured = [data[k] for k in measured_names]\n",
    "\n",
    "# Switch terms\n",
    "gamma_f,gamma_r = data['gamma_f'],data['gamma_r']\n",
    "\n",
    "# DUT\n",
    "dut_meas = data['DUT']\n",
    "\n",
    "# 50 ohm termination\n",
    "res_50ohm_meas = data['res_50ohm']"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple Multiline TRL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the line lengths in meters (including thru)\n",
    "l = [0, 0.3e-3, 2.3e-3]\n",
    "\n",
    "# Do the calibration\n",
    "cal = NISTMultilineTRL(\n",
    "    measured = measured,  # Measured standards\n",
    "    Grefls = [-1], # Reflection coefficient of the reflect, -1 for short\n",
    "    l = l,         # Lengths of the lines\n",
    "    er_est = 7,    # Estimate of transmission line effective permittivity\n",
    "    switch_terms = (gamma_f, gamma_r) # Switch terms\n",
    "    )\n",
    "\n",
    "# Correct the DUT using the above calibration\n",
    "corrected = cal.apply_cal(dut_meas)\n",
    "\n",
    "corrected.plot_s_db()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare calibrations with different combinations of lines"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we loop through different line combinations to demonstrate the difference in the calibration accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run NIST Multiline TRL calibration with different combinations of lines\n",
    "\n",
    "# Put through and reflect to their own list ...\n",
    "mtr = measured[:2]\n",
    "\n",
    "# and lines on their own\n",
    "mlines = measured[2:]\n",
    "\n",
    "line_len = l[1:]\n",
    "\n",
    "cals = []\n",
    "duts = []\n",
    "\n",
    "line_combinations = [[0], [1], [0,1]]\n",
    "\n",
    "for used_lines in line_combinations:\n",
    "\n",
    "    m = mtr + [mlines[i] for i in used_lines]\n",
    "\n",
    "    # Add thru length to list of line lengths\n",
    "    l = [l[0]] + [line_len[i] for i in used_lines]\n",
    "\n",
    "    # Do the calibration\n",
    "    cal = NISTMultilineTRL(\n",
    "        measured = m,  # Measured standards\n",
    "        Grefls = [-1], # Reflection coefficient of the reflect, -1 for short\n",
    "        l = l,         # Lengths of the lines\n",
    "        er_est = 7,    # Estimate of transmission line effective permittivity\n",
    "        switch_terms = (gamma_f, gamma_r) # Switch terms\n",
    "        )\n",
    "\n",
    "    # Correct the DUT using the above calibration\n",
    "    corrected = cal.apply_cal(dut_meas)\n",
    "    corrected.name = f'DUT, lines {used_lines}'\n",
    "\n",
    "    duts.append(corrected)\n",
    "    cals.append(cal)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transmission of the corrected DUT\n",
    "Plot the corrected DUT calibrated with different combination of calibration lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.title('DUT S21')\n",
    "for dut in duts:\n",
    "    dut.plot_s_db(m=1, n=0)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### S11 of corrected DUT with different amount of calibration lines\n",
    "\n",
    "S11 shows bigger changes. \n",
    "\n",
    "* With one short line low frequencies are very noisy\n",
    "* With only the long line the calibration is very inaccurate at frequencies where the phase difference of the thru and line is close to a multiple of 180 degrees\n",
    "* With both lines calibration accuracy is good everywhere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.title('DUT S11')\n",
    "for dut in duts:\n",
    "    dut.plot_s_db(m=0, n=0)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalized standard deviation of different calibrations\n",
    "\n",
    "Normalized standard deviation can be used to measure the accuracy of the calibration. Lower number means calibration is less sensitive to the measurement noise.\n",
    "\n",
    " * TRL calibration with one 90 degrees long line has normalized standard deviation of 1. \n",
    " * TRL calibration with one 180 degree long lossless line is singular and has infinite normalized standard deviation.\n",
    " * With multiple lines normalized standard deviation less than one is possible.\n",
    " \n",
    "Note that the nstd is normalized such that it doesn't consider the actual measurement noise. It's calculated only from the solved propagation constant and line lengths. The threshold of how large it can be depends on the DUT being measured, measurement noise and the required accuracy of the measurement. If there are large spikes, such as are visible in the long line case below, that's a sign that the calibration is very close to singular at that frequency and the measurement accuracy is going to be poor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_ghz = dut.frequency.f_scaled\n",
    "\n",
    "plt.figure()\n",
    "plt.title('Calibration normalized standard deviation')\n",
    "for e, cal in enumerate(cals):\n",
    "    plt.plot(f_ghz, cal.nstd, label=f'Lines: {line_combinations[e]}')\n",
    "plt.ylim([0,20])\n",
    "plt.legend(loc='upper right')\n",
    "dut.frequency.labelXAxis()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate effective complex relative permittivity of transmission lines used in the calibration\n",
    "\n",
    "Effective complex relative permittivity $\\epsilon_{r,eff}$ of a transmission line is related to the propagation constant $\\gamma$ as:\n",
    "\n",
    "$\\gamma = \\frac{2\\pi f}{c}\\sqrt{\\epsilon_{r,eff}}$, where $c$ equals the speed of light and $f$ is frequency.\n",
    "\n",
    "In general it's a complex value with the imaginary part indicating losses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define calibration standard media\n",
    "freq = dut.frequency\n",
    "\n",
    "# Get the cal with the both lines\n",
    "cal = cals[-1]\n",
    "\n",
    "plt.figure()\n",
    "plt.title('CPW effective permittivity (real part)')\n",
    "plt.plot(f_ghz, cal.er_eff.real, label='Solved er_eff')\n",
    "plt.xlabel('Frequency (GHz)')\n",
    "plt.legend(loc='lower right')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TRL calibration accuracy is the best when line length difference is 90 degrees. Solved propagation constant and effective permittivity however are more accurate the bigger the line length difference is. At low frequencies the estimate is noisier due to the line phase difference being small."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the phase of the solved reflection coefficient"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying the calibration to the measured reflect standard we can get the calibrated S-parameters of the unknown reflect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.title('Solved reflection coefficient of the reflect standard')\n",
    "cal.apply_cal(measured[1]).plot_s_deg(n=0, m=0, label='Solved short')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reference plane shift\n",
    "\n",
    "Because propagation constant of the media is solved during the calibration it's possible to shift the reference plane by a specified distance.\n",
    "\n",
    "The reference plane shift can be specified with `ref_plane` argument. The shift should be specified in meters, negative lengths is towards the VNA. By default the same shift is applied to both ports. Unequal shift on the two ports is supported by passing a two element list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cal_shift = NISTMultilineTRL(\n",
    "    measured = measured,  # Measured standards\n",
    "    Grefls = [-1], # Reflection coefficient of the reflect, -1 for short\n",
    "    l = l,         # Lengths of the lines\n",
    "    er_est = 7,    # Estimate of transmission line effective permittivity\n",
    "    switch_terms = (gamma_f, gamma_r), # Switch terms\n",
    "    # Shift reference planes twords VNA by this amount (in m) on both ports\n",
    "    ref_plane = -50e-6\n",
    "    )\n",
    "\n",
    "# Correct the DUT using the above calibration\n",
    "corrected_thru = cal.apply_cal(measured[0])\n",
    "corrected_thru_shifted = cal_shift.apply_cal(measured[0])\n",
    "\n",
    "corrected_thru.plot_s_deg(m=1, n=0, label='Thru phase')\n",
    "corrected_thru_shifted.plot_s_deg(m=1, n=0, label='Reference plane shifted thru phase')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calibration reference impedance renormalization\n",
    "\n",
    "The reference impedance of the calibration is by default the transmission line characteristic impedance. If we know the actual characteristic impedance of the lines we can give it to the calibration routine with the `z0_line` argument to renormalize the measured S-parameters to a fixed reference `z0_ref`.\n",
    "\n",
    "If the conductance per unit length (G) is much lower than the capacitive reactance per unit length ($j\\omega C_0$), the characteristic impedance of the transmission line can be written in terms of the propagation constant $\\gamma$ and capacitance per unit length $C_0$:\n",
    "\n",
    "$Z_0 = \\gamma/(j 2 \\pi f C_0)$\n",
    "\n",
    "If $C_0$ is known it can be given to the calibration routine with `c0` parameter to renormalize the calibration reference impedance to `z0_ref` (defaults to 50 ohms) assuming G = 0.\n",
    "If the line is lossy the characteristic impedance is complex valued and giving a single `c0` instead of a fixed `z0_line` is usually more accurate.\n",
    "\n",
    "In this case we know that the line characteristic impedance is actually 55 ohms. To renormalize the calibration from 55 ohms to 50 ohms we need to give `z0_line=55` argument to the calibration routine. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cal_ref = NISTMultilineTRL(\n",
    "    measured = measured,  # Measured standards\n",
    "    Grefls = [-1], # Reflection coefficient of the reflect, -1 for short\n",
    "    l = l,         # Lengths of the lines\n",
    "    er_est = 7,    # Estimate of transmission line effective permittivity\n",
    "    switch_terms = (gamma_f, gamma_r), # Switch terms\n",
    "    z0_line = 55, # Line actual characteristic impedance\n",
    "    z0_ref = 50 # Calibration reference impedance\n",
    "    )\n",
    "\n",
    "cal.apply_cal(res_50ohm_meas).s11.plot_s_db(label=r'50 $\\Omega$ termination |$S_{11}$|, Z_ref = line')\n",
    "cal_ref.apply_cal(res_50ohm_meas).s11.plot_s_db(label=r'50 $\\Omega$ termination |$S_{11}$|, Z_ref = 50 $\\Omega$')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After renormalization the 50 ohm termination measurement shows good matching. It's not perfectly matched due to the noise in the measurements."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TUG multiline TRL\n",
    "\n",
    "TUG multiline TRL (`TUGMultilineTRL`) is a method for implementing multiline TRL calibration using a different approach than that of NIST. The definition of the calibration standards in TUG multiline is the same as in NIST multiline, requiring at least two lines and symmetric reflect. More mathematical details about TUG multiline can be found here: https://ziadhatab.github.io/posts/multiline-trl-calibration/\n",
    "\n",
    "Below is an example comparing NIST multiline and TUG multiline. Both methods deliver the same answer. However, later in the [Monte Carlo (MC) analysis](#monte-carlo-analysis---nist-vs-tug), we will see that the TUG method delivers better statistical performance compared to the NIST method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load all measurement data into a dictionary\n",
    "data = skrf.io.general.read_all_networks('multiline_trl_data/')\n",
    "# Pull out measurements by name into an ordered list\n",
    "measured_names = ['thru','reflect','linep3mm','line2p3mm']\n",
    "measured = [data[k] for k in measured_names]\n",
    "# Switch terms\n",
    "gamma_f,gamma_r = data['gamma_f'],data['gamma_r']\n",
    "# DUT\n",
    "dut_meas = data['DUT']\n",
    "# define the line lengths in meters (including thru)\n",
    "l = [0, 0.3e-3, 2.3e-3]\n",
    "\n",
    "# NIST multiline TRL\n",
    "cal = NISTMultilineTRL(\n",
    "    measured = measured,  # Measured standards\n",
    "    Grefls = [-1], # Reflection coefficient of the reflect, -1 for short\n",
    "    l = l,         # Lengths of the lines\n",
    "    er_est = 7,    # Estimate of transmission line effective permittivity\n",
    "    switch_terms = (gamma_f, gamma_r) # Switch terms\n",
    "    )\n",
    "\n",
    "# Correct the DUT using the above calibration\n",
    "corrected = cal.apply_cal(dut_meas)\n",
    "corrected.name = 'DUT NIST'\n",
    "\n",
    "corrected.plot_s_db()\n",
    "\n",
    "# TUG multiline TRL - note the difference in\n",
    "cal = TUGMultilineTRL(\n",
    "        line_meas=[measured[0]]+measured[2:],\n",
    "        line_lengths=l,\n",
    "        er_est=7,\n",
    "        reflect_meas=[measured[1]],\n",
    "        reflect_est=[-1],\n",
    "        switch_terms = (gamma_f, gamma_r)\n",
    "        )\n",
    "\n",
    "# Correct the DUT using the above calibration\n",
    "corrected = cal.apply_cal(dut_meas)\n",
    "corrected.name = 'DUT TUG'\n",
    "\n",
    "corrected.plot_s_db(linestyle='--')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calibration without reflect measurement\n",
    "\n",
    "One nice feature of `TUGMultilineTRL` is that you can use only line standards if you only care about the propagation constant and S21 and S12 of the calibrated DUT. However, to obtain correct S11 and S22 values of the calibrated DUT, you also need to provide measurement of the reflect standard. This can be seen in the example below, where S11 and S22 do not overlap when using and not using the reflect standard."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Full TUG multiline TRL\n",
    "cal = TUGMultilineTRL(\n",
    "        line_meas=[measured[0]]+measured[2:],\n",
    "        line_lengths=l,\n",
    "        er_est=7,\n",
    "        reflect_meas=[measured[1]],\n",
    "        reflect_est=[-1],\n",
    "        switch_terms = (gamma_f, gamma_r)\n",
    "        )\n",
    "\n",
    "# Correct the DUT using the above calibration\n",
    "corrected = cal.apply_cal(dut_meas)\n",
    "corrected.name = 'TUG w/ reflect'\n",
    "\n",
    "corrected.plot_s_db(linestyle='-')\n",
    "\n",
    "# TUG multiline without reflect measurement\n",
    "cal = TUGMultilineTRL(\n",
    "        line_meas=[measured[0]]+measured[2:],\n",
    "        line_lengths=l,\n",
    "        er_est=7,\n",
    "        switch_terms = (gamma_f, gamma_r)\n",
    "        )\n",
    "\n",
    "# Correct the DUT using the above calibration\n",
    "corrected = cal.apply_cal(dut_meas)\n",
    "corrected.name = 'TUG w/o reflect'\n",
    "\n",
    "corrected.plot_s_db(linestyle='--')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calibration eigenvalue\n",
    "\n",
    "The NIST method incorporates a normalized standard deviation to emphasize the sensitivity of the calibration. Similarly, the TUG method employs a metric based on the eigenvalue of the calibration problem. The TUG multiline method solves only one weighted eigendecomposition, resulting in a real-valued eigenvalue $\\lambda$. This eigenvalue signifies the stability of the calibration. The closer the eigenvalue is to zero, the more sensitive the calibration is. In essence, it can be viewed as the opposite of the normalized standard deviation defined in the NIST method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_ghz = dut.frequency.f_scaled\n",
    "\n",
    "plt.figure()\n",
    "plt.title('TUG multiline eigenvalue')\n",
    "plt.plot(f_ghz, cal.lambd)\n",
    "dut.frequency.labelXAxis()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monte Carlo analysis - NIST vs TUG\n",
    "\n",
    "Using MC analysis, we generate synthetic data of the standards with random noise and repeat the process multiple times, each time with a new noise sample. At each trial, we also compute and store the calibrated DUT, which we use to calculate the mean absolute error (MAE) of the calibration method. The general definition of MAE is as follows:\n",
    "\n",
    "$$\n",
    "\\text{MAE} = \\frac{1}{N}\\sum_{i=1}^{N} |x_i - x_\\mathrm{true}|\n",
    "$$\n",
    "\n",
    "where $x_i$ is the $i$-th MC trial and $x_\\mathrm{true}$ is the true value. The code below generates the reference true value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "freq = skrf.Frequency(1, 100, 201, 'GHz')\n",
    "# CPW media used for DUT and the calibration standards\n",
    "cpw = CPW(freq, w=40e-6, s=51e-6, ep_r=12.9, t=5e-6, rho=2e-8)\n",
    "# 1.0 mm coaxial media for calibration error boxes\n",
    "coax1mm = Coaxial(freq, z0_port=50, Dint=0.434e-3, Dout=1.0e-3, sigma=1e8)\n",
    "f_ghz = cpw.frequency.f*1e-9\n",
    "# Error boxes\n",
    "X = coax1mm.line(1, 'm', z0=58, name='X')\n",
    "Y = coax1mm.line(1.1, 'm', z0=40, name='Y')\n",
    "# Realistic looking switch terms\n",
    "gamma_f = coax1mm.delay_load(0.2, 21e-3, 'm', z0=60)\n",
    "gamma_r = coax1mm.delay_load(0.25, 16e-3, 'm', z0=56)\n",
    "# Lengths of the lines used in the calibration, units are in meters\n",
    "line_lengths = [0, 0.3e-3, 2.3e-3]  # first is thru\n",
    "lines = [cpw.line(l, 'm') for l in line_lengths]\n",
    "# Attenuator with mismatched feed lines\n",
    "dut_feed = cpw.line(1.5e-3, 'm', z0=60)\n",
    "dut = dut_feed**cpw.attenuator(-10)**dut_feed\n",
    "# Reflect standard\n",
    "reflect_offset = 10e-6\n",
    "short = cpw.delay_short(reflect_offset, 'm')\n",
    "short = two_port_reflect(short, short)\n",
    "\n",
    "# Measured... embedded in the error boxes\n",
    "line_meas = [terminate(X**k**Y,gamma_f, gamma_r) for k in lines]\n",
    "reflect_meas = [terminate(X**short**Y, gamma_f, gamma_r)]\n",
    "dut_meas = terminate(X**dut**Y, gamma_f, gamma_r)\n",
    "er_est = 6.3-0.0001j\n",
    "reflect_est = [-1]\n",
    "\n",
    "# NIST multiline TRL (Noiseless case as reference)\n",
    "measured = [line_meas[0]] + [reflect_meas[0]] + line_meas[1:]\n",
    "cal = NISTMultilineTRL(\n",
    "    measured = measured,\n",
    "    Grefls = reflect_est,\n",
    "    l = line_lengths,\n",
    "    refl_offset = reflect_offset,\n",
    "    er_est = er_est,\n",
    "    switch_terms=[gamma_f, gamma_r])\n",
    "\n",
    "dut_ideal = cal.apply_cal(dut_meas)\n",
    "dut_ideal.name = 'DUT'\n",
    "dut_ideal.plot_s_db()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now start the MC analysis by introducing noise in polar coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Monte Carlo Analysis\n",
    "print('MC starts: ')\n",
    "M = 100 # number of trials\n",
    "mag_std = 0.05 # std of magnitude (linear)\n",
    "phase_std = 20 # std of phase (degrees)\n",
    "\n",
    "cal_dut_NIST = []\n",
    "cal_dut_TUG  = []\n",
    "\n",
    "for inx in range(M):\n",
    "    # add noise (make sure to copy and not override the original data)\n",
    "    lines_n   = [NW.copy() for NW in line_meas]\n",
    "    for NW in lines_n:\n",
    "        NW.add_noise_polar(mag_std, phase_std)\n",
    "\n",
    "    reflect_n   = [NW.copy() for NW in reflect_meas]\n",
    "    for NW in reflect_n:\n",
    "        NW.add_noise_polar(mag_std, phase_std)\n",
    "\n",
    "    # TUG multiline TRL\n",
    "    cal = TUGMultilineTRL(line_meas=lines_n, line_lengths=line_lengths, er_est=er_est,\n",
    "                reflect_meas=reflect_n, reflect_est=reflect_est, reflect_offset=reflect_offset,\n",
    "                switch_terms=[gamma_f, gamma_r])\n",
    "    cal_dut_TUG.append(cal.apply_cal(dut_meas))\n",
    "\n",
    "    # NIST multiline TRL\n",
    "    measured = [lines_n[0]] + [reflect_n[0]] + lines_n[1:]\n",
    "    cal = NISTMultilineTRL(\n",
    "        measured = measured,\n",
    "        Grefls = reflect_est,\n",
    "        l = line_lengths,\n",
    "        refl_offset = reflect_offset,\n",
    "        er_est = er_est,\n",
    "        switch_terms=[gamma_f, gamma_r])\n",
    "    cal_dut_NIST.append(cal.apply_cal(dut_meas))\n",
    "    print(f'MC Index: {inx+1}/{M} done.', end='\\r', flush=True)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we compute and plot the MAE. The results show that TUG multiline has a lower error compared to the NIST method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dut_MAE(MC, ideal):\n",
    "    # compute mean absolute error\n",
    "    return np.array([ abs(x.s-ideal.s) for x in MC ]).mean(axis=0)\n",
    "\n",
    "MAE_NIST = dut_MAE(cal_dut_NIST, dut_ideal)\n",
    "MAE_TUG  = dut_MAE(cal_dut_TUG, dut_ideal)\n",
    "\n",
    "f_ghz = dut_ideal.frequency.f_scaled\n",
    "\n",
    "plt.figure()\n",
    "plt.title('Mean Absolute Error')\n",
    "plt.plot(f_ghz, MAE_NIST[:,0,0], label='NIST, S11')\n",
    "plt.plot(f_ghz, MAE_NIST[:,1,0], label='NIST, S21')\n",
    "plt.plot(f_ghz, MAE_TUG[:,0,0], label='TUG, S11')\n",
    "plt.plot(f_ghz, MAE_TUG[:,1,0], label='TUG, S21')\n",
    "plt.legend(loc='upper right')\n",
    "dut_ideal.frequency.labelXAxis()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation to generate the input data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is how we made the data used above. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create frequency and Media "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "freq = skrf.Frequency(1, 100, 201, 'GHz')\n",
    "\n",
    "# CPW media used for DUT and the calibration standards\n",
    "cpw = CPW(freq, w=40e-6, s=51e-6, ep_r=12.9,\n",
    "                     t=5e-6, rho=2e-8)\n",
    "print(cpw.z0[0])\n",
    "\n",
    "# 1.0 mm coaxial media for calibration error boxes\n",
    "coax1mm = Coaxial(freq, z0_port=50, Dint=0.434e-3, Dout=1.0e-3, sigma=1e8)\n",
    "print(coax1mm.z0[0])\n",
    "\n",
    "f_ghz = cpw.frequency.f*1e-9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make realistic looking error networks.\n",
    "\n",
    "Propagation constant determination is iterative and doesn't work as well when the error networks are randomly generated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = coax1mm.line(1, 'm', z0=58, name='X')\n",
    "Y = coax1mm.line(1.1, 'm', z0=40, name='Y')\n",
    "\n",
    "plt.figure()\n",
    "plt.title('Error networks')\n",
    "X.plot_s_db()\n",
    "Y.plot_s_db()\n",
    "\n",
    "# Realistic looking switch terms\n",
    "gamma_f = coax1mm.delay_load(0.2, 21e-3, 'm', z0=60)\n",
    "gamma_r = coax1mm.delay_load(0.25, 16e-3, 'm', z0=56)\n",
    "\n",
    "plt.figure()\n",
    "plt.title('Switch terms')\n",
    "gamma_f.plot_s_db()\n",
    "gamma_r.plot_s_db()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Fictitious measurements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lengths of the lines used in the calibration, units are in meters\n",
    "line_len = [0.3e-3, 2.3e-3]\n",
    "lines = [cpw.line(l, 'm') for l in line_len]\n",
    "\n",
    "# Attenuator with mismatched feed lines\n",
    "dut_feed = cpw.line(1.5e-3, 'm', z0=60)\n",
    "dut = dut_feed**cpw.attenuator(-10)**dut_feed\n",
    "\n",
    "res_50ohm = cpw.resistor(50) ** cpw.short(nports=2) ** cpw.resistor(50)\n",
    "\n",
    "# Through and non-ideal short\n",
    "# Real reflection coefficient is solved during the calibration\n",
    "\n",
    "short = cpw.delay_short(10e-6, 'm')\n",
    "\n",
    "actuals = [\n",
    "    cpw.thru(),\n",
    "    two_port_reflect(short, short),\n",
    "    ]\n",
    "\n",
    "actuals.extend(lines)\n",
    "\n",
    "# Measured\n",
    "measured = [X**k**Y for k in actuals]\n",
    "\n",
    "# Switch termination\n",
    "measured = [terminate(m, gamma_f, gamma_r) for m in measured]\n",
    "\n",
    "# Add little noise to the measurements\n",
    "for m in measured:\n",
    "    m.add_noise_polar(0.001, 0.1)\n",
    "\n",
    "names = ['thru', 'reflect', 'linep3mm', 'line2p3mm']\n",
    "for k,name in enumerate(names):\n",
    "    measured[k].name=name\n",
    "\n",
    "\n",
    "# Noiseless DUT so that all the noise will be from the calibration\n",
    "dut_meas = terminate(X**dut**Y, gamma_f, gamma_r)\n",
    "dut_meas.name = 'DUT'\n",
    "\n",
    "res_50ohm_meas = terminate(X**res_50ohm**Y, gamma_f, gamma_r)\n",
    "res_50ohm_meas.name = 'res_50ohm'\n",
    "\n",
    "# Put through and reflect to their own list ...\n",
    "mtr = measured[:2]\n",
    "\n",
    "# and lines on their own\n",
    "mlines = measured[2:]\n",
    "\n",
    "# write data to disk\n",
    "write_data = False\n",
    "if write_data:\n",
    "    [k.write_touchstone(dir='multiline_trl_data/') for k in measured]\n",
    "    gamma_f.write_touchstone('multiline_trl_data/gamma_f.s1p')\n",
    "    gamma_r.write_touchstone('multiline_trl_data/gamma_r.s1p')\n",
    "    dut_meas.write_touchstone(dir='multiline_trl_data/')\n",
    "    res_50ohm_meas.write_touchstone(dir='multiline_trl_data/')"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}